#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <gmp.h>
#include <time.h>
#include<sys/time.h>
#define INDEX(i,j,k) (i+j)%k
#define INDEX2(i,j) INDEX(i,j,2)
#define INDEX3(i,j) INDEX(i,j,3)
#define NUMDIG 10000001


/* Método de Gauss-Legendre sequêncial */
int main()
{
    mpf_t precisao;  
    mpf_t pi;
    mpf_t a[2], b[2], t[2], p[2];
    mpf_t old_pi;
    int i, i0, i1;
    //VARIAVEIS PARA GUARDAR O TEMPO
    int tempo;
    struct timeval inicio, fim;
       
    gmp_printf("calculando precisao... \n");
    /* Precisão em bits = ceil( NDIGITS * log2(10) ) */
    mpf_set_default_prec(33219281);
    //mpf_set_default_prec(33223);
    mpf_init_set_ui(precisao, 10);
    mpf_pow_ui(precisao, precisao, NUMDIG); 
    mpf_ui_div(precisao, 1, precisao); // precisao = 1/(10^PREC)
    //gmp_printf("precisao = %F\n\n", precisao);
    mpf_init(pi);    

    gmp_printf("inicializando valores... \n");    
    
    /* PEGA O TEMPO INICIAL DE EXECUCAO */
    gettimeofday(&inicio, NULL);
    
    /* Calcula a0 = 1 */
    mpf_init_set_ui(a[0], 1);

    /* Calcula b0 = 1/sqrt(2) */
    mpf_init_set_ui(b[0], 2); mpf_sqrt(b[0], b[0]); mpf_ui_div(b[0], 1, b[0]);

    /* Calcula t0 = 1/4 */
    mpf_init_set_d(t[0], 0.25);

    /* Calcula p0 = 1 */
    mpf_init_set_ui(p[0], 1);

    mpf_inits(a[1], b[1], t[1], p[1], old_pi, NULL);
    gmp_printf("calculando pi... \n");    
    for(i = 0; ; i++)
    {
      //gmp_printf("Iteracao %i\n", i);
      i0 = INDEX2(i, 0);
      i1 = INDEX2(i, 1);
      
      /* Calcula os parâmetros */
      mpf_add(a[i1], a[i0], b[i0]); 
      mpf_div_ui(a[i1], a[i1], 2); // a1 = (a0+b0)/2
      mpf_mul(b[i1], a[i0], b[i0]); 
      mpf_sqrt(b[i1], b[i1]); // b1 = sqrt(a0*b0)
      mpf_sub(t[i1], a[i0], a[i1]); 
      mpf_mul(t[i1], t[i1], t[i1]); 
      mpf_mul(t[i1], t[i1], p[i0]); 
      mpf_sub(t[i1], t[i0], t[i1]); // t1 = t0-p0*(a0-a1)^2
      mpf_mul_ui(p[i1], p[i0], 2); // p1 = 2*p0
      mpf_set(old_pi, pi);
      //gmp_printf("old_pi[%i] => %Ff\n", i, old_pi);
      
      /* Calcula o valor de pi = (a1+b1)^2/(4*t1)*/
      mpf_add(pi, a[i1], b[i1]); 
      mpf_mul(pi, pi, pi); 
      mpf_div(pi, pi, t[i1]); 
      mpf_div_ui(pi, pi, 4); 
      //gmp_printf("pi[%i] = %Ff\n", i, pi);
      
      /* Verifica convergência */
      mpf_t aux;
      int convergencia;
      mpf_init(aux);

      /* aux = |pi-old_pi| */
      mpf_sub(aux, pi, old_pi); 
      //gmp_printf("aux = %Ff\n", aux);
      mpf_abs(aux, aux);
      
      /* Compara com a precisão e, se for menor, garantimos convergência */	
      convergencia = (mpf_cmp(aux, precisao) <= 0);
      //gmp_printf("convergencia = %i\n\n", convergencia);
      mpf_clear(aux); 
      if(convergencia) break;
    }
    
    mpf_clears(a[0], a[1], b[0], b[1], t[0], t[1], p[0], p[1], old_pi, NULL);
 
    gmp_printf("numero de iteracoes: %i \nvalor estimado de pi: %.1000Ff \n",i,pi);
    
    //PEGA O TEMPO FINAL DE EXECUCAO
    gettimeofday(&fim, NULL);

    //DIFERENCA ENTRE FINAL E INICIAL
    tempo = fim.tv_sec - inicio.tv_sec;
    printf("TEMPO DE EXECUÇÃO %d segundos\n", tempo);

    gmp_printf("imprimindo arquivo 'pi-gauss-legendre-seq.txt'\n");
    FILE *fp = fopen("pi-gauss-legendre-seq.txt", "w"); /* Abre o arquivo */
    mpf_out_str(fp, 10, NUMDIG, pi); /* Imprime com precisão 10^7 dígitos*/   
    gmp_printf("\n");
    fclose(fp); 
    
    mpf_clear(precisao);
    mpf_clear(pi);
    return 0;
}
